import os
import sys

from astropy.io import fits as pyfits
from astropy.visualization.interval import ZScaleInterval
from skimage import data, exposure, img_as_float
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from astropy.visualization import simple_norm
# import Image,ImageFont,ImageDraw
import skimage as image
from pylab import *
import numpy as np
from astropy.wcs import WCS
import PIL.Image as pil


def coordinate_transform(datafile, xmin, xmax, ymin, ymax, newdatafile):
    xff = open(newdatafile, 'w')
    #print(f'xshift={xmin}')
    #print(f'yshift={ymin}')
    for line in open(datafile,'r'):
        values = [s for s in line.split()]
        if float(values[0]) > xmin and float(values[0]) < xmax:
            if float(values[1]) > ymin and float(values[1]) < ymax:
                xc = float(values[0]) - xmin - 1.0
                yc = float(values[1]) - ymin  - 1.0
                xff.write("%f %f\n" % (yc,xc))

    xff.close()


def circle_points(RN, center, radius):
    radians = np.linspace(0, 2 * np.pi, RN)
    #print("============")
    #print(center[1])
    #print(center[0])
    c = center[1] + radius * np.cos(radians)  # polar co-ordinates
    r = center[0] + radius * np.sin(radians)
    #print(np.array([c, r]).T)
    return np.array([c, r]).T

def plotstar(filename, markcolor, radius):
    markcolor = markcolor
    datafile = filename
    if os.path.exists(datafile):
        ff = open(datafile, 'r')
        newdataall = ff.readlines()
        #print("newdataall=%s" % (newdataall))
        #print("To read the para for the plot")
        for f1 in newdataall:
            f1t = f1.split()
            # print(f1t)
            xs = float(f1t[0])
            ys = float(f1t[1])
            # print(xs)
            # print(ys)
            # radius = 10
            points = circle_points(200, [xs, ys], radius)[:-1]
            plt.plot(points[:, 0], points[:, 1], color=markcolor, linestyle="solid", lw=0.3)


def calRotZero(filfits):
    # !/usr/bin/env python
    # - to cal rotation angle for a fits image.
    # - example: calRotZero.py [fitsfile_wcscalibrated]
    # --
    from mpdaf.obj import WCS
    from astropy.io import fits
    import sys
    import os

    hdulist = fits.open(filfits)
    w = WCS(hdulist[0].header)
    rot = w.get_rot()
    # print (f"{filfits} rot_angle = {rot} deg ")
    print(f"{filfits}  {rot} ")
    # print (f"{rot}" )
    return rot

def cutImageA(fullimage, xmin, xmax, ymin, ymax, otname, xc, yc, error_arcmin, comment, instrument,
              devname, rot, filter,ra,dec,exptime, limitmag):
    newimg_data = []
    fullimage = fullimage
    print(fullimage)
    xmin = xmin
    xmax = xmax
    ymin = ymin
    ymax = ymax
    otname = otname
    xc = xc
    yc = yc
    error_arcmin = error_arcmin
    comment = comment
    instrument = instrument
    devname = devname
    rot =  rot
    filter = filter
    limitmag = "%.2f" % float(limitmag)
    exptime = exptime
    ra = ra
    dec = dec
    pixelscale = 0.7   #default value
    if instrument == "CGFT":
        pixelscale = 0.62  #3C-CCD-g, 3C-CCD-r, 3C-CCD-i
        if instrument == "CGFT" and devname == "PF-CMOS15":
            pixelscale = 1.4 #arcseconds
        if instrument == "CGFT" and devname == "PF-CMOS12":
            pixelscale = 1.2 #arcseconds, needed to varify
    if instrument == "GWAC-F60A" or instrument == "GWAC-F60B":
        pixelscale = 0.52
    #pixelscale = 0.62 #arcsec for CGFT, and 0.56 for F60
    radius = error_arcmin * 60 / pixelscale # pixelscale
    print(pixelscale)
    print(radius)
    #print(fullimage)
    #print(xmin)
    #print(xmax)
    #print(ymin)
    #print(ymax)

    # cutimage = "aa.png"
    #cutimage = fullimage[0:-4] + "_" + otname + ".png"
    cutimage = fullimage[0:-4]  + ".png"
    cutimage_norotate = fullimage[0:-4]  + "_norotate.png"
    cutimagefit = fullimage[0:-4] + "_" + otname + ".fit"
    print(f' cutimage is {cutimage}')

    with pyfits.open(fullimage) as hdul:
        #hdul.info()
        hdul.verify('fix')
        newimage = hdul[0]
        img_data = np.array(hdul[0].data)
        #print(img_data.shape)
        wcs = WCS(hdul[0].header)
        #print(wcs)
        #print(img_data)
        print(xmin)
        print(xmax)
        print(ymin)
        print(ymax)
        # print(len(img_data))
        # print(img_data.shape)
        #newimg_data = img_data[xmin:xmax, ymin:ymax]
        newimg_data = img_data[ymin:ymax, xmin:xmax]
        #print(newimg_data.shape)
        # print(newimg_data)
       # print("@@@@@@@@@@@@@@@@@")
       # print(len(newimg_data))
       # print(newimg_data.shape)

        pyfits.writeto(cutimage, newimg_data, header=None, overwrite=True)

        #================================================
        #to get the subimage in fit format
        #newimage.data = img_data[xmin:xmax, ymin:ymax]
        newimage.data = img_data[ymin:ymax, xmin:xmax]
        #print(newimage.shape)

        newimage.header = hdul[0].header
        phdr = newimage.header
        CCDSEC="[%d:%d,%d:%d]" % (xmin, xmax, ymin,ymax)
        BIASSEC="[%d:%d,%d:%d]" % (1,xmax-xmin+1, 1, ymax-ymin+1)
        TRIMSEC=CCDSEC
        print(TRIMSEC)
        WCSDIM = 2
        LTM1_1 = 1
        LTM2_2 = 1
        WAT0_001='system=image'
        WAT1_001='wtype=tan axtype=ra'
        WAT2_001='wtype=tan axtype=dec'
        LTV1=-(xmin-1)
        LTV2=-(ymin-1)
        #print(LTV1)
        #print(LTV2)
        CRPIX1=phdr['CRPIX1']+LTV1
        CRPIX2=phdr['CRPIX2']+LTV2
        #print( CRPIX1)
        #print(CRPIX2)
        phdr.set('WCSDIM',WCSDIM)
        phdr.set('LTM1_1',LTM1_1)
        phdr.set('LTM2_2',LTM2_2)
        phdr.set('WAT0_001',WAT0_001)
        phdr.set('WAT1_001',WAT1_001)
        phdr.set('WAT2_001',WAT2_001)
        phdr.set('CRPIX1',CRPIX1)
        phdr.set('CRPIX2',CRPIX2)
        phdr.set('CCDSEC',CCDSEC)
        phdr.set('BIASSEC',BIASSEC)
        phdr.set('LTV1',LTV1)
        phdr.set('LTV2',LTV2)
        #phdr.set('TRIMSEC',TRIMSEC)
        #phdr.set('IMWHOLE',fullimage)
        phdr.set('EXTEND',False)
        #print(phdr)
        newimage.writeto(cutimagefit,overwrite=True)
        #============================

        interval = ZScaleInterval()
        vmin, vmax = interval.get_limits(newimg_data)
        print(f"vmin={vmin}")
        print(f'vmax={vmax}')
        #points = circle_points(200, [xc, yc], radius)[:-1]
        xbox = (xmax - xmin) / 2
        ybox = (ymax - ymin) / 2
        print(f'xbox={xbox} and ybox={ybox}')

        points = circle_points(200, [xbox, ybox], radius)[:-1]
        # plt.subplot(1, 1, 1)
        #fig = plt.gcf()
        #plt.axis('off')
        #fig.set_size_inches(7.0 / 3, 7.0 / 3)
        #plt.gca().xaxis.set_major_locator(plt.NullLocator())
        #plt.gca().yaxis.set_major_locator(plt.NullLocator())
        #plt.subplots_adjust(top=1, bottom=0, left=0, right=1, hspace=0, wspace=0)
        #plt.margins(0, 0)
        print("to show the wcs")
        #plt.subplot(projection=wcs)
        #ax = plt.gca()
        #ax.set_aspect(1)
        #plt.figure()
        plt.grid(True, color='white', linestyle='--', linewidth=0.3)
        print("to show the imshow")
        #vmin = vmin
        vmax = vmax - 20
        #print(newimg_data)
        plt.imshow(newimg_data, cmap=plt.cm.gray, vmin=vmin, vmax=vmax, origin='lower', aspect='equal')

        #plt.subplots_adjust(top=1, bottom=0, right=1, left=0, hspace=0, wspace=0)
        #plt.margins(0, 0)

        datafile = "PS1_xy.txt"
        if os.access(datafile,os.F_OK):
            newdatafile = "PS1_xy_shift.txt"
            coordinate_transform(datafile, xmin, xmax, ymin, ymax, newdatafile)
            plotstar(filename=newdatafile, markcolor="yellow", radius=5)

        datafile = "newOT_xy.txt"
        if os.access(datafile,os.F_OK):
            newdatafile = "newOT_xy_shift.txt"
            coordinate_transform(datafile, xmin, xmax, ymin, ymax, newdatafile)
            plotstar(filename=newdatafile, markcolor="red", radius=5)
            plotstar(filename=newdatafile, markcolor="red", radius=20)

        datafile = "usno_xy.txt"
        if os.access(datafile,os.F_OK):
            newdatafile = "usno_xy_shift.txt"
            coordinate_transform(datafile, xmin, xmax, ymin, ymax, newdatafile)
            plotstar(filename=newdatafile, markcolor="cyan", radius=7)
        # datafile = "MXT.txt"
        # if os.access(datafile,os.F_OK):
        #     newdatafile = "MXT_shift.txt"
        #     coordinate_transform(datafile, xmin, xmax, ymin, ymax, newdatafile)
        #     plotstar(filename=datafile, markcolor="cyan", radius=468)


        # plt.axis('off')
        # plt.savefig(fname=cutimage_norotate, format='png', dpi=200, bbox_inches='tight', pad_inches=0.0)
        # #plt.subplot()
        # print("to show the imshow")
        # #im1 = plt.imshow(newimg_data, cmap=plt.cm.gray, vmin=vmin, vmax=vmax, origin='lower', aspect='equal')
        # #rot = "80"   #定义了旋转角度，目前CGFT主焦点是旋转量，先上下镜像，再转转
        # #rot = calRotZero(fullimage)
        # im1 = pil.open(cutimage_norotate)
        # #im1 = im1.transpose(pil.FLIP_TOP_BOTTOM)
        # #im3 = im2.rotate(360 - eval(str(rot)))
        # im1 = im1.rotate(0 - eval(str(rot)))
        # plt.imshow(im1)
        plt.axis('on')
        print("to show the circle")
        plt.plot(points[:, 0], points[:, 1], color='white', linestyle="dashed", lw=1)
        plt.title(f"r={error_arcmin}arcmin, {comment}", fontsize=10)
        textmsg1 = str(filter) + " " + "limitmag=" + str(limitmag)
        textmsg2 = "R.A.=" + str(ra) + ",Dec.=" + str(dec)
        textmsg3 = "exptime=" + str(exptime) + "sec"
        textmsg4 = textmsg2 + "," + textmsg3 + "," + textmsg1
        #plt.text(10, 10, textmsg, )
        plt.text(10, 10, textmsg4, fontsize=7, color="blue",
                 bbox=dict(facecolor='white', alpha=1))
        # plt.text(25, 10, textmsg2, fontsize=7, color="blue",
        #          bbox=dict(facecolor='white', alpha=1))
        # plt.text(10, 10, textmsg1, fontsize=7, color="blue",
        #          bbox=dict(facecolor='white', alpha=1))
        plt.savefig(fname=cutimage, format='png', dpi=500, bbox_inches='tight', pad_inches=0.0)
        #plt.show()

        #plt.show()


if __name__ == "__main__":
    print(f'length of args is {len(sys.argv)}')
    if len(sys.argv) == 18:
        print('aaa')
        # img = "G171013_C05906_0696.fit"
        img = str(sys.argv[1])
        xmin = str(sys.argv[2])
        xmax = str(sys.argv[3])
        ymin = str(sys.argv[4])
        ymax = str(sys.argv[5])
        triggername = str(sys.argv[6])
        xc = str(sys.argv[7])
        yc = str(sys.argv[8])
        errorarcmin = float(sys.argv[9])
        obstime = str(sys.argv[10])
        instrument = str(sys.argv[11])
        obsfilter = str(sys.argv[12])
        devname = str(sys.argv[13])  # 终端名称，PF-CMOS15, PF-CMOS12, ...
        rot = str(sys.argv[14])
        exptime = str(sys.argv[15])
        ra = str(sys.argv[16])
        dec = str(sys.argv[17])
        limitmag = "99.9"
        # img = "data/sb22113056_g_1670072664_0001_1.fit"
        # xmin = 1378.9
        # xmax = 1564.6
        # ymin = 2525.3
        # ymax = 2711.0
        # triggername = "sb22"
        # xc = 1471.8
        # yc = 2618.2
        # errorarcmin = 1.0
        # obstime = "2022-12-03T13:04:50.1"
        # instrument = "CGFT"
        # obsfilter = "g"
        # devname = "PF-CMOS15"
        comment = instrument + " " + obsfilter + "_" + obstime
        print(img)
        print(xmin)
        print(xmax)
        print(ymin)
        print(ymax)
        print(xc)
        print(yc)
        print(errorarcmin)
        print(instrument)
        print(devname)
        print(rot)
        print(limitmag)
        print(comment)
        cutImageA(fullimage=img, xmin=int(float(xmin)), xmax=int(float(xmax)),
                  ymin=int(float(ymin)), ymax=int(float(ymax)),
                  otname=triggername,
                  xc=int(float(xc)), yc=int(float(yc)),
                  error_arcmin=errorarcmin, comment=comment,
                  instrument=instrument, devname=devname, rot=float(rot), filter=obsfilter, limitmag=limitmag,
                  ra=ra,dec=dec,exptime=exptime)

    else:
        print('bbb')
        # img = "G171013_C05906_0696.fit"
        img = str(sys.argv[1])
        xmin = str(sys.argv[2])
        xmax = str(sys.argv[3])
        ymin = str(sys.argv[4])
        ymax = str(sys.argv[5])
        triggername = str(sys.argv[6])
        xc = str(sys.argv[7])
        yc = str(sys.argv[8])
        errorarcmin = float(sys.argv[9])
        obstime = str(sys.argv[10])
        instrument = str(sys.argv[11])
        obsfilter = str(sys.argv[12])
        devname = str(sys.argv[13])  # 终端名称，PF-CMOS15, PF-CMOS12, ...
        rot = str(sys.argv[14])
        exptime = str(sys.argv[15])
        ra = str(sys.argv[16])
        dec = str(sys.argv[17])
        limitmag = str(sys.argv[18])
        # img = "data/sb22113056_g_1670072664_0001_1.fit"
        # xmin = 1378.9
        # xmax = 1564.6
        # ymin = 2525.3
        # ymax = 2711.0
        # triggername = "sb22"
        # xc = 1471.8
        # yc = 2618.2
        # errorarcmin = 1.0
        # obstime = "2022-12-03T13:04:50.1"
        # instrument = "CGFT"
        # obsfilter = "g"
        # devname = "PF-CMOS15"
        comment = instrument + " " + obsfilter + "_" + obstime
        print(img)
        print(xmin)
        print(xmax)
        print(ymin)
        print(ymax)
        print(xc)
        print(yc)
        print(errorarcmin)
        print(instrument)
        print(devname)
        print(rot)
        print(limitmag)
        print(comment)
        cutImageA(fullimage=img, xmin=int(float(xmin)), xmax=int(float(xmax)),
                  ymin=int(float(ymin)), ymax=int(float(ymax)),
                  otname=triggername,
                  xc=int(float(xc)), yc=int(float(yc)),
                  error_arcmin=errorarcmin, comment=comment,
                  instrument=instrument, devname=devname, rot=float(rot), filter=obsfilter, limitmag=limitmag,
                  ra=ra,dec=dec,exptime=exptime)



